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ABSTRACT 

Hill's equations are an approximation that is useful in a number of areas of 
astrophysics including planetary rings and planetesimal disks. We derive a sym- 
plectic method for integrating Hill's equations based on a generalized leapfrog. 
This method is implemented in the parallel iV-body code, PKDGRAV and tested 
on some simple orbits. The method demonstrates a lack of secular changes in 
orbital elements, making it a very useful technique for integrating Hill's equa- 
tions over many dynamical times. Furthermore, the method allows for efficient 
collision searching using linear extrapolation of particle positions. 



Subject headings: methods: iV-body simulations - methods: numerical 
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Introduction 



There are a number of situations in planetary dynamics that require the exploration 
of near-circular orbits . Cur rent topics of interest i n this category inc l ude planetary rings 



( I Wisdom &: Tremaind Il988l ) and planet formation (ITanga et al.l 120041 ; iBarnes et al.l 12009 



hereafter BQLR). In this circumsta nce t he equ ations of motion can be linearized about the 
circular orbit as was first done by iHilll (118781 ) to study the lunar orbit. That is, the mo- 
tions of bodies are described with respect to a Cartesian frame that is in uniform circular 
motion about a central body, and excursions from the center of the frame are small com- 
pared to the distance to the central body. In the absence of perturbations, the resulting 
Hill's equations describe simple epicyclic motiom and can also be used for disk dynamics 
( IGoldreich &l Lvnden-Belllll965l ; lJulian fc Toomrdll966l ) and the escape of stars from globu- 
lar clusters (jHeggidl200ll ). 



A noticeable disadvantage of Hill's equations for numerical integration is that they 
contain a velocity-dependent force . Simulation codes for large iV-body simulations, (e.g. 
Springell 120051 ; IWadsley et al.ll2004j ) typically use the leapfrog integration scheme, which is 
second order, symplectic, and easy to implement. The leapfrog scheme can be modified 
to take velocity-dependent forces into account and still retain second order, as is done for 
Smoothed Particle Hydrodynamics (SPH); however, this destroys its symplectic nature. 

The power of symplectic integrators is rooted in the property that any truncation error 
can be represented as a perturbing Hamiltonian. Hence for sufficiently small step size, the 
numerical system has conserved quantities similar to the integrals of motion of the physical 
system. That is, the numerical integration is an exact solution to an approximate Hamilto- 
nian. This property is particularly important when following systems for many dynamical 
times such as the long-term evolution of the Solar System, or investigating the stability 
of extrasolar planetary systems. In these situations, if the integrator introduces secular 
changes in the actions, the dynamics being investigated can be fund amentally changed. 
Hence symplectic integrators are widely used in such investiga t ions (IHolman fc Wisdom 
19931 : iLevison fc Duncan! 1 199^ iMalhotral fl~995l : iLee fc Peald I2OO2I ; lllivera fc Lissauerl boOph T 
Both planetary rings and planetesimal dynamics are systems that evolve over large numbers 
of dynamical times, and therefore may also benefit from the use of symplectic integrators. 



A symplectic integrator for Hill's equations was introduced by iHeggid (1200 ll ) in the 
context of escape of stars from globular clusters. However, in that work, the integrator is 
not actually put to use; instead, the orbits were calculated using a Hermite integrator. The 
symplectic integrator was expressed as an implicit set of equations which, however, could be 
solved explicitly. As shown below, a canonical transformation can significantly simplify the 
Hamiltonian, and therefore simplify the resulting integrator. 
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Saha fc Tremaind (119921 ) introduced a formalism (also see IWisdom fc Holmanl Il99ll ) 



for deriving symplectic integrators of a generalized leapfrog type by separating the Hamil- 
tonian into parts that can be integrated exactly and then using com mutator algebra to 
combine these solutions into a symplectic solution to the full problem. iQuinn et al.l (119971 ) 
showed how this technique could be used for cosmological simulations that involve a time- 
depende nt Hamiltonian. This formalism has also been used to construct higher-order inte- 
grators ( Laskar fc Robutel I2OOI ; Chambers fc Murison 2000 ). integrators that handle close 



encounters ( Chambers 



1999: 



Duncan et al. 



1998h . and integrators that efficiently in tegrate 



problems with a large dynamic range (jSaha &: Tremaindll994l ; iMcNeil fc Nelsonll2009i ). Here 
we will apply the technique to Hill's equations. In section 2 present the Hamiltonian formula- 
tion of Hill's equations from which in section 3 we derive a symplectic integrator suitable for 
use in a large TV-body code . In section 4, we describe its implementation in the PKDGRAV 
TV-body code (jStadell l200ll ) . explicitly stating the algorithm for timestepping a simulation, 
and in section 5, we perform tests appropriate for the application of planetesimal dynamics 
in the early Solar System. Section 6 contains a short discussion and summary. 



Hamiltonian Formulation 



The Lagrangian for Hill's equations in the orbital plane is 



C=±[(x- Sly) 2 + (y + Qx) 2 ] + l -tf(2x 2 - y 2 ) - *(x, y), (1) 



( jHeggid 120011 ) where x and y are, respectively, the distances perpendicular to and along the 
direction of rotation from the center of a frame in circular motion with angular speed Q. $ 
is the potential due to other forces, e.g., interactions with other particles. In all that follows, 
we will neglect the motion in the z direction since it is trivial to integrate in the standard 
way. 

Lagrange's equations give the standard Hill's equations of motion, 

x-2n y -3n 2 x = (2) 

y + 2Qx = -—. (3) 
dy 

In this form, the presence of the velocity-dependent terms requires a modification to the 
leapfrog method such that a predicted velocity is used in the estimate of the final acceleration. 
This maintains second order, but is obviously not time reversible and destroys the symplectic 
nature of leapfrog. However, as we show below, a symplectic integrator can be derived for 
this system. 
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To construct a symplectic integrator, we first derive the Hamiltonian form of the equa- 
tions of motion. From their definitions, the canonical momenta are 

Px = =x-ily (4) 
Py = -qT = V + fix (5) 



and the Hamiltonian is 



H(x,y,p x ,p y ) = ^ + P l + n(yp x - xp y ) - ifi 2 (2x 2 - y 2 ) + <S>{x,y). (6) 

Now consider a new set of canonical coordinates, (X, Y, P x , P y ), derived from the generating 
function 

S 2 (x,y,P x ,P y ) = xP x + yP y -ttxy. (7) 
The rules of canonical transformations then give 

Px = ^ = Px- ty; p y = = P y - fix, (8) 

and 

In terms of the original positions and velocities, these new canonical coordinates are x, y, 
P x = x, and Py = y + 2Qx. The Hamiltonian in these coordinates is 

P 2 P 2 Q 2 x 2 
H(x, y, P x , P y ) = ^ + ^-- 2QxP y + — + y), (10) 

a somewhat simpler form than equation (151) . 

Hamilton's equations of motion are therefore 

x = P x (11) 
y = P y -2ttx (12) 

P x = 2QP y - fi 2 x - ^ (13) 

From these equations, it is obvious that P y is constant in the absence of perturbing forces. 
This is equivalent to the conservation of angular momentum, and leads to the conserved 
quantity, Yli^yi m a many-particle system. In particular, Y^i^y * s conserved in a collision 
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between particles ( IWisdom fc Tremaindll988l ). which will be useful when calculating collision 
outcomes (see below). However, for periodic boundary conditions, P y will change as a particle 
crosses the boundary in x because of the shear across the box. Nevertheless the time-averaged 
total P v should be constant for a system that does not have a net motion in the x direction 
( IWisdom fc Tremaindll988l ). Also from the equations of motion in this form it is clear why 
Hill's equations are easy to integrate numerically. If the guiding center of the motion is at 
x = 0, then P y = for all time in the absence of perturbations, and the motion reduces to 
a harmonic oscillator with frequency Q. 



3. Symplectic Integrators 

Normally one can cr eate a symplectic i ntegr ator by separating the Hamiltonian into 



exactly integral parts as in lSaha fc Tremainei (119921 ). but the presence of a velocity-dependent 



force makes this approach nontrivial in the case of Hill's problem. The Hamiltonian can be 
split as follows, 

H = \{P 2 X + P 2 y ) + -2nxP y + ^- + Hx, y), (15) 



where Hpp is the free particle Hamiltonian and Hmdf is considered to be the (momentum- 
dependent) "force" term. The Hpp Hamiltonian is easily solved and is just the motion of a 
particle with constant velocity: 

x(t Q + r) = x(t )+rP x (t ) 

y(t Q + r) = y(t )+rP y (t ) (16) 
P x (t + r) = P x (t ) 

Py{t + r) = Py{t ), 

where r is the timestep and to is the initial time. However, the equation of motion corre- 
sponding to the Hmdf Hamiltonian can not be solved easily. This is because for this part of 
the Hamiltonian y is not constant (OHmdf I dP y ^ 0), so one must evaluate the force along 
a trajectory determined by y = —2Qx to solve for P y , and then use this to solve for P x . This 
would prove intractable in a large simulation. 

Instead, let us separate the Hamiltonian into a mixed term, Hmix, and a momentum- 
independent force term, Hmif, as follows: 

H = \{P 2 X + Py 2 ) - 2QxP y + ^ + $(*, y) . (17) 



2 v x y, y 2 

v ' V v 

Hmix Hmif 
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Hmif is easily integrated to give the equations of motion, 



P x (t + t) 
P y (t + r) 



<9$ 

p x (t ) - t ( n 2 x(t ) + — 



to, 



Py(to) 



<9$ 

dy 



to 



The mixed Hamiltonian, Hmixi gives the equations of motion, 

x P x 



y 

A. 



P y - 2Qx 

2QPy 

0. 



(19) 
(20) 
(21) 
(22) 



These can be integrated exactly as follows. P y is a constant so P y {t) = P y (t ). P x (t) can 
now be solved. x(t) can be solved once P x (t) is known, and finally y(t) can be solved since 
we know P y (t) and x(t). We therefore have 



x(t + r) = x(t ) + rP x (t ) + r 2 nP y (t ) 

y(to + r) = y(t ) + r(Py(t )-2Qx(t ))~T 2 QP x (t )-r 3 ^n 2 P y (t ) 

P,(to + r) = P x (t ) + r2nP y (t ) 

Py(t + r) = P y (t ). 



(23) 



These can be used to construct a second-order symplectic integrator exactly analo- 
gous to leapfrog by applying equations (TT51) for half a timestep, equations f[2"B"j) for a full 
timestep, and equations (fl8l) for another half timestep. If we were simply integrating force 
equations this is straightforward to implement in a large iV-body code. However, in the 
case of planetesimal and planetary ring dynamics, collisions between particles need to be 
detec ted. Current collision d etection algorithms rely on the position updates being linear in 



time ([Richardson et al.l 120001 ). and certainly not cubic in time as in the above. 



In an attempt to simplify the mixed equations of motion, let us separate Huix even 
further into the free particle Hamiltonian, Hfp, and the cross term, He = —2QxP y . That 
is, 

H = \{P 2 X + P 2 y ) + ~2QxP y + ?^ + «&(*, y) . (24) 




Hfp Hmif 
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The cross Hamiltonian is easily integrable giving 

x(t + r) = x(t ) (25) 

y(t + r) = y(t ) - r2VLx(t ) (26) 

P x (t + r) = P x (t ) + T2{lP y (t ) (27) 

P y (t + r) = P y (t ). (28) 

Therefore, a first-order symplectic scheme presents itself as follows. 1) Update the momenta 
using equations (ITS]) . 2) Update P x using equation (I2"7j) and the P y from step 1. 3) Update 
the y positions of the particles according to equation (125]) . If collisions are being considered, 
they are searched for in this step. 4) Perform a standard position update using the free 
particle Hamiltonian (equation ITS"]) , again searching for collisions. In our implementation, 
steps 3 and 4 are combined into a single position update that includes the collision search. 



The construction of a second-order scheme follows using the formalism of lSaha &: Tremaine 



(]1992l ). If we refer to the evolution of phase space for a time r under the Hamiltonians Hmif, 
H FP , and H c as MIF{r) (eq. [T8j), FP{t) (eq. Mj, and C(r) (eq. I25H28J respectively, then 
the combination of operators MIF(t/2)C(t/2)FP(t)C(t/2)MIF(t/2) will evolve the sys- 
tem for a timestep r with second-order acc uracy. That is, the error Hamiltonian will be of 
order r 3 or higher (jSaha fc Tremaine! Il992l ). 



4. Implementation in an iV-body code 



To test the usefulness of this formulation we have implem ented the above integration 
algorithm in the parallel gravity code PKDGRAV (jStadell 120011 ) as part of the technique for 



solving an iV-body system in a patch corotating in a Kepler potential ( iRichardson et al. 2000 



Porco et al. 20081 ; iBQLPj ). This code uses a standard form of leapfrog, where the velocities 
are first updated by a half step, a Kick, then the positions are updated by a full step, a Drift, 
and finally the velocites are given a second half step Kick. Only minor changes were needed 
to implement the above second-order algorithm. The most straightforward way to modify 
the algorithm is to change the Kick routine so that it only includes the terms present in 
the Hmif Hamiltonian, and modify the Drift routine to include the operations of equations 
(125]) through (128]) as well as the standard Drift of the positions. However, the existing Drift 
routine in PKDGRAV is complicated by the handling of periodic boundary conditions and 
the search for collisions, so we instead rearranged the operations so that the Drift remains a 
simple linear extrapolation of the positions with constant velocities. 

In detail, the modifications are as follows. After the opening Kick routine updates 
the velocities according to Hmif, h calculates the canonical momentum, P y , and updates x 
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(which is equivalent to P x ) using equation fl27|) . However, y is updated to be the sum of 
the contributions of the free particle Hamiltonian, Hpp, the cross term (1261) applied at the 
beginning of the Drift (for time r/2), and the cross term applied at the end of the Drift (also 
for a time r/2): 



(Also see equation [331 below.) This y along with x can now be used to linearly update 
the positions according to both "Cross" operators and the free particle operator using an 
essentially unmodified Drift routine. Finally, the closing Kick again uses equation ff27|) to 
update x, sets y to be P y — 2Qx, and updates the velocities according to Hmif- The above 
algorithm requires storage for a new attribute, P y , for each particle. 

One modification to the Drift routine involves the handling of periodic boundary con- 
ditions. Often, Hill's equations are integrated with periodic boundaries in the x and y 
directions. If a particle exits the computational volume in the, e.g. , positive y direction, 
then it is replaced by a particle with the same velocity and x coordinate on the negative 
y boundary. Handling the x boundaries is a little more complicated because of the shear 
across the patch: a particle in a circular orbit on the outer boundary is moving slower in the 
y direction than a particle on a circular orbit on the inner boundary by an amount |fiAx, 
where Ax is the width of the patch. This implies an increase in the P y of the particle of 
^QAx. This corresponds to the fact that in a Kepler potential, the circular velocity of an or- 
bit decreases outwards, while the angular momentum of a circular orbit increases outwards. 
Hence when the Drift routine detects a particle has crossed the x boundary, P y is changed 
accordingly. 

The Drift routine also performs collision detection and resolution. Therefore if a par- 
ticle's momentum is changed due to a collision, or if a new particle is created either due to 
merging or fragmentation, the canonical momentum, P y , needs to be updated to reflect the 
change. Solving equation (1291) for P y provides a simple means to calculate a new P y at any 
time during the drift. For a merger, the conservation of the total P y of the bodies involved 
in the collision could be used to assign a P y to the merged particle. No other changes to 
PKDGRAV were needed to implement this algorithm. 

In summary, a single timestep of a single particle starting from the state (x n , y n ,i n , y n ) 




(29) 
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and producing the state (x n+ i, y n +x, x n +\, Vn+i) is performed as follows. 

Xn+l/A = X n — - \il X n ~\ — 1 (30) 

P, = fc + Hfa, (31) 

2 ay 

in+i/2 = x n+1/i + rQP y (32) 

Vn+l/2 = Py- ~ Q(x n + TX n+1 / 2 ) (33) 

X n+1 = X n + TX n+ i/2 (34) 

Vn+x = yn + ry n+1 /2 (35) 

in+3/4 = in+1/2 + TQPy (36) 

x n+1 = x n+3/4 - - I il x n+ i H — 1 (37) 

p on Td$(x n+1 ,y n+1 ) 

y n+1 = P y -2ilx n+ i-- — (38) 

Although fractional subscripts are used for intermediate values of the state variables for 
clarity, the only extra storage needed for this update is for the canonical momentum, P y . 

Equations ([30] - [33]) are implemented in the first Kick routine, and equations ([36] - [38]) are 
implemented in the last Kick routine, leaving a simple form for the Drift, equations (I34p and 
( 1351) . during which boundary crossing and collision detection is performed. 



5. Tests of the Method 

As a first test of the usefulness of a symplectic scheme for Hill's equations, we used the 
implementation in PKDGRAVto integrate a single particle in a Kepler potential and compare 
the conserved integrals of the system using our new integrator with those using a standard 
second-order integration method. The standard method integrates equations ([2]) and ([3]) 
with the Kick- Drift- Kick leapfrog described above, with the modification that velocities are 
predicted to the end of the timestep using the old accelerations in order to calculate the 
velocity dependent part of the force. This is the same algorithm that is used to handle th e 



velocity dependent forces arising in Smooth Particle Hydrodynamics (jWadsley et al.ll2004[ ). 



In this case of a single particle $ = 0, and a physically relevant combination of the 
integrals is the eccentricity, which in terms of the canonical coordinates can be expressed as 

e2 = p^ " mxPy + fi V + Ap y"> ' (39) 
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Fig. 1. — Relative change in eccentricity for a single particle. Ae/e is plotted vs orbit 
number for a 100-orbit integration of a particle with an initial eccentricity of 0.001. The 
dotted line shows the evolution using the standard non-symplectic integration algorithm 
with 100 steps per orbit. The solid line shows the evolution using the symplectic algorithm 
with 20 steps per orbit. The symplectic algorithm with 100 steps per orbit has a maximum 
Ae/e of 0.00017, so its evolution would be a horizontal line at the resolution of this plot. All 
integrations were performed with PKDGRAV on a single processor. 
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where R is the radius of the orbit of the patch. Figure [T] shows the relative change in 
eccentricity in an integration of a particle with an initial e = 0.001 for 100 orbits. The 
dotted line shows this change for the standard integration method when 100 steps per orbit 
are used, while the solid line shows results for our new method with only 20 steps per orbit. 
The new integrator demonstrates the typical behavior of a symplectic method: the numerical 
value of the integral of motion oscillates around the true value, and there is no secular drift. 
An integration with the symplectic integrator using 100 steps per orbit has a maximum Ae/e 
of 0.00017. Comparing this with the maximum Ae/e for the integration with 20 steps per 
orbit, 0.0043, indicates that the error is scaling as r 2 , where r is the timestep. This is as 
expected for a second-order integrator. A careful inspection of the figure shows that the 
dotted line has a slope that is slightly increasing with time, implying that in the standard 
method, the eccentricity drift in this case grows faster for larger eccentricities. 

A somewhat more relevant test of the integrator is following an encounter in the re- 
stricted three body approximation. Specifically, we use the PKDGRAV implementation to 
integrate the orbit of a test particle as it comes within a Hill radius of a massive body on a 
circular orbit. Such a situation is not uncommon in simulations of planetesimal growth: the 
large bodies are in somewhat circular orbits (BQLR). The accuracy of the integration can 



be evaluated using the Jacobi integral (IDuncan et al 



19891 ) 



r = 3QV-, 2 -^ + (x2 2 + ^ )1/2 , (40) 

where m is the mass of the massive body. Guided by the end-state configuration of the L\ 
simulation in BQLR, we set the mass of the massive body to be 3.78 x 10 18 g (300 times 
the mass of a 1 km planetesimal). The test particle is placed in an orbit such that it 
comes within one Hill radius of the massive body at aphelion, and the encounter speed at 
closest approach is given by the RMS velocity in the BQLR simulation, 2 m s _1 . These 
parameters imply an eccentricity of 1.6 x 10~ 4 and a relative difference in semi-major axis 
e = \a — a m \/a = 1.69 x 10 -4 , where a m is the semi-major axis of the massive body. We 
follow the motion of the test particle starting at perihelion, through the conjunction with 
the massive body and to the subsequent perihelion. Due to the encounter, the eccentricity of 
the test body changes by 6.4 x 10~ 7 . This is somewhat gre ater than the eccentricity change 



expected from the mapping formula of IDuncan et al.l (119891 ) . 1.5 x 10 7 , presumably because 



this encounter does not satisfy their approximation that e ^ e. 

Figure [2] shows how well T is conserved during this encounter as a function of the 
integration timestep for the symplectic and the non-symplectic integrations. As with the 
eccentricity in the simple orbit case, the error in the Jacobi constant scales as r 2 . However, 
for a given timestep the symplectic integration algorithm gives an order of magnitude im- 
provement in the conservation of T. At the largest timestep plotted, r = P/30, the test 
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r/P 

Fig. 2. — Maximum relative change in the Jacobi constant, T, during an encounter with a 
massive body as a function of step size. A single orbit during which the test body comes 
within one Hill radius of a massive body is integrated between two successive perihelia. 
The dotted line shows the maximum change in T over this orbit as a function of step size 
in units of the orbital period for the standard non-symplectic integration. The solid line 
shows the same quantity for the symplectic algorithm. All integrations were performed with 
PKDGRAV on a single processor. 
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body moves through the encounter at about 2 ru per timestep. Hence reasonably resolving 
the encounter requires r = P/100 or smaller. 



6. Discussion and Summary 

The choice of algorithm for a numerical simulation is critical to obtaining accurate phys- 
ical results. This is particularly true for simulations performed over many dynamical times 
where there is an opportunity for truncation error to build up in a secular manner. In this 
case, an algorithm that appears to work well over a few orbits may produce unacceptably 
incorrect results when used for hundreds of orbits. Using the standard second-order integra- 
tor to solve Hill's equations in the context of planetesimal accretion illustrates this problem. 
Although individual orbits are followed reasonably well with a few hundred steps per orbit, 
the secular growth in eccentricity shown in Fig. [1] could overwhelm any physical changes 
in the eccentricity distribution of planetesimals in a hundred orbits or so. Faced with this 
problem, one must either go to a higher-order algorithm, which is difficult to implement in 
a large parallel simulation code, or use much smaller timesteps, which significantly increases 
the computational expense. 

Fortunately, for the case of Hill's equation we have discovered a second-order sym- 
plectic integrator that does not display any secular growth in eccentricity. Moreover, our 
solution is linear in the time extrapolation of particle positions, permitting efficient collision 
detection, which is valuable for most large N simulat ions of planetesim a ls and planetary 



rings. The algorithm is derived from the formalism of ISaha fc Tremaind (119921 ) where the 
Hamiltonian is separated into parts that by themselves are integrable. This separation is in 
turn made possible via a canonical transformation to coordinates that significantly simpli- 
fies the Hamiltonian. The algorithm has been implemented in the scalable parallel iV-body 
code PKDGRAV, and this code is currently being used for follow-on simulations to those 
described in BQLR, which simulated the growth of planetesimals over hundreds of orbits. 

As described above, our integrator has a fixed timestep which is inefficient for simulations 
that have to resolve close encounters between bodies where the encounter timescales are a 
small fraction of the orbital time. For example, in planet formation simulations the encounter 
timescale is hours compared to an orbital time of one year. Constructing an integrator that 
can adjust ti mesteps in order to handle close encoun ters and yet retains symplectic properties 



requires care (IDuncan et al.lll998l ; IChamberslll999l ). In our implementation, we have simply 



adjusted the timestep of particles that are experiencing a collision. Although this destroys 
the symplectic properties of the integrator, for the planetesimal simulations in BQLR, a 
typical particle makes about one hundred orbits before experiencing a collision. Hence, 
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even with the non-symplectic timestep adjustment, our integrator significantly improves the 
quality of these simulations. 
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